Fit spatiotemporal diet models with covariates
Max Lindmark
05 December, 2022
Intro
In this script, I take the collated stomach data set and calculate aggregates (feeding ratio, total weight of prey groups) and predictor variables for diet data, aggregate to get 1 stomach = 1 row per prey type (not prey individual). I also select only the columns I need for model fitting, join environmental covariates and cpue covariates for cod and flounder, and lastly saduria biomass densities.
Load packages & source functions
# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(readxl)
library(tidylog)
library(RCurl)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(patchwork)
library(viridis)
library(sdmTMB)
library(mapplots)
library(sp)
library(raster)
library(ggeffects)
library(modelr)
# Source code for lan_lot to UTM
source("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/R/functions/lon_lat_utm.R")
# Source code for map plots
source("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/R/functions/map_plot.R")
# Trim the map plot to better fit stomach data
xmin2 <- 303000
xmax2 <- 916000
xrange <- xmax2 - xmin2
ymin2 <- 5980000
ymax2 <- 6580000
yrange <- ymax2 - ymin2
theme_set(theme_plot())
plot_map_fc2 <- plot_map_fc +
theme(legend.position = "bottom") +
xlim(xmin2*1.5, xmax2*0.9) +
ylim(ymin2*1.025, ymax2*0.98268)
# Continuous colors
options(ggplot2.continuous.colour = "viridis")Read stomach data
small_cod_stomach <- read_csv("data/clean/small_cod_stomach.csv") %>%
drop_na(oxy, temp, depth, scod_kg_km2, lcod_kg_km2, fle_kg_km2) %>%
rename(density_cod_small = scod_kg_km2,
density_cod_large = lcod_kg_km2,
density_flounder = fle_kg_km2) %>%
mutate(fyear = as.factor(year),
fquarter = as.factor(quarter),
ices_rect = as.factor(ices_rect),
oxy_sc = (oxy - mean(oxy)) / sd(oxy),
temp_sc = (temp - mean(temp)) / sd(temp),
depth_sc = (depth - mean(depth)) / sd(depth),
density_saduria_sc = (density_saduria - mean(density_saduria)) / sd(density_saduria),
density_cod_small_sc = (density_cod_small - mean(density_cod_small)) / sd(density_cod_small),
density_cod_large_sc = (density_cod_large - mean(density_cod_large)) / sd(density_cod_large),
density_flounder_sc = (density_flounder - mean(density_flounder)) / sd(density_flounder),
haul_id = as.factor(haul_id),
species2 = "Small cod")
#> Rows: 1153 Columns: 33
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (8): pred_id, predator_latin_name, species, ices_rect, haul_id, pred_we...
#> dbl (25): pred_weight_g, pred_length_cm, year, quarter, month, day, subdiv, ...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> drop_na: no rows removed
#>
#> rename: renamed 3 variables (density_flounder, density_cod_large, density_cod_small)
#>
#> mutate: converted 'ices_rect' from character to factor (0 new NA)
#>
#> converted 'haul_id' from character to factor (0 new NA)
#>
#> new variable 'fyear' (factor) with 6 unique values and 0% NA
#>
#> new variable 'fquarter' (factor) with 2 unique values and 0% NA
#>
#> new variable 'oxy_sc' (double) with 152 unique values and 0% NA
#>
#> new variable 'temp_sc' (double) with 152 unique values and 0% NA
#>
#> new variable 'depth_sc' (double) with 64 unique values and 0% NA
#>
#> new variable 'density_saduria_sc' (double) with 74 unique values and 0% NA
#>
#> new variable 'density_cod_small_sc' (double) with 158 unique values and 0% NA
#>
#> new variable 'density_cod_large_sc' (double) with 157 unique values and 0% NA
#>
#> new variable 'density_flounder_sc' (double) with 159 unique values and 0% NA
#>
#> new variable 'species2' (character) with one unique value and 0% NA
large_cod_stomach <- read_csv("data/clean/large_cod_stomach.csv") %>%
drop_na(oxy, temp, depth, scod_kg_km2, lcod_kg_km2, fle_kg_km2) %>%
rename(density_cod_small = scod_kg_km2,
density_cod_large = lcod_kg_km2,
density_flounder = fle_kg_km2) %>%
mutate(fyear = as.factor(year),
fquarter = as.factor(quarter),
ices_rect = as.factor(ices_rect),
oxy_sc = (oxy - mean(oxy)) / sd(oxy),
temp_sc = (temp - mean(temp)) / sd(temp),
depth_sc = (depth - mean(depth)) / sd(depth),
density_saduria_sc = (density_saduria - mean(density_saduria)) / sd(density_saduria),
density_cod_small_sc = (density_cod_small - mean(density_cod_small)) / sd(density_cod_small),
density_cod_large_sc = (density_cod_large - mean(density_cod_large)) / sd(density_cod_large),
density_flounder_sc = (density_flounder - mean(density_flounder)) / sd(density_flounder),
haul_id = as.factor(haul_id),
species2 = "Large cod")
#> Rows: 2153 Columns: 33
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (8): pred_id, predator_latin_name, species, ices_rect, haul_id, pred_we...
#> dbl (25): pred_weight_g, pred_length_cm, year, quarter, month, day, subdiv, ...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> drop_na: no rows removed
#>
#> rename: renamed 3 variables (density_flounder, density_cod_large, density_cod_small)
#>
#> mutate: converted 'ices_rect' from character to factor (0 new NA)
#>
#> converted 'haul_id' from character to factor (0 new NA)
#>
#> new variable 'fyear' (factor) with 6 unique values and 0% NA
#>
#> new variable 'fquarter' (factor) with 2 unique values and 0% NA
#>
#> new variable 'oxy_sc' (double) with 164 unique values and 0% NA
#>
#> new variable 'temp_sc' (double) with 164 unique values and 0% NA
#>
#> new variable 'depth_sc' (double) with 67 unique values and 0% NA
#>
#> new variable 'density_saduria_sc' (double) with 76 unique values and 0% NA
#>
#> new variable 'density_cod_small_sc' (double) with 158 unique values and 0% NA
#>
#> new variable 'density_cod_large_sc' (double) with 171 unique values and 0% NA
#>
#> new variable 'density_flounder_sc' (double) with 171 unique values and 0% NA
#>
#> new variable 'species2' (character) with one unique value and 0% NA
flounder_stomach <- read_csv("data/clean/flounder_stomach.csv") %>%
drop_na(oxy, temp, depth, scod_kg_km2, lcod_kg_km2, fle_kg_km2) %>%
rename(density_cod_small = scod_kg_km2,
density_cod_large = lcod_kg_km2,
density_flounder = fle_kg_km2) %>%
mutate(fyear = as.factor(year),
fquarter = as.factor(quarter),
ices_rect = as.factor(ices_rect),
oxy_sc = (oxy - mean(oxy)) / sd(oxy),
temp_sc = (temp - mean(temp)) / sd(temp),
depth_sc = (depth - mean(depth)) / sd(depth),
density_saduria_sc = (density_saduria - mean(density_saduria)) / sd(density_saduria),
density_cod_small_sc = (density_cod_small - mean(density_cod_small)) / sd(density_cod_small),
density_cod_large_sc = (density_cod_large - mean(density_cod_large)) / sd(density_cod_large),
density_flounder_sc = (density_flounder - mean(density_flounder)) / sd(density_flounder),
haul_id = as.factor(haul_id),
species2 = "Flounder")
#> Rows: 2579 Columns: 33
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (8): pred_id, predator_latin_name, species, ices_rect, haul_id, pred_we...
#> dbl (25): pred_weight_g, pred_length_cm, year, quarter, month, day, subdiv, ...
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> drop_na: no rows removed
#>
#> rename: renamed 3 variables (density_flounder, density_cod_large, density_cod_small)
#>
#> mutate: converted 'ices_rect' from character to factor (0 new NA)
#>
#> converted 'haul_id' from character to factor (0 new NA)
#>
#> new variable 'fyear' (factor) with 6 unique values and 0% NA
#>
#> new variable 'fquarter' (factor) with 2 unique values and 0% NA
#>
#> new variable 'oxy_sc' (double) with 153 unique values and 0% NA
#>
#> new variable 'temp_sc' (double) with 153 unique values and 0% NA
#>
#> new variable 'depth_sc' (double) with 64 unique values and 0% NA
#>
#> new variable 'density_saduria_sc' (double) with 77 unique values and 0% NA
#>
#> new variable 'density_cod_small_sc' (double) with 148 unique values and 0% NA
#>
#> new variable 'density_cod_large_sc' (double) with 160 unique values and 0% NA
#>
#> new variable 'density_flounder_sc' (double) with 161 unique values and 0% NA
#>
#> new variable 'species2' (character) with one unique value and 0% NA
# Load cache
# qwraps2::lazyload_cache_dir(path = "R/main_analysis/main_diet_analysis_cache/html")Read the prediction grid (scales variables wrt data!)
# I need predicted cod and flounder densities for this
# pred_grid_1 <- read_csv("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/data/clean/pred_grid_(1_2).csv")
# pred_grid_2 <- read_csv("/Users/maxlindmark/Dropbox/Max work/R/cod_interactions/data/clean/pred_grid_(2_2).csv")
#
# pred_grid <- bind_rows(pred_grid_1, pred_grid_2)Plot data in space
full_dat <- bind_rows(small_cod_stomach, large_cod_stomach, flounder_stomach) %>%
pivot_longer(c(tot_feeding_ratio, saduria_feeding_ratio, benthic_feeding_ratio),
names_to = "prey_group",
values_to = "feeding_ratio")
#> pivot_longer: reorganized (tot_feeding_ratio, benthic_feeding_ratio, saduria_feeding_ratio) into (prey_group, feeding_ratio) [was 5885x43, now 17655x42]
## Total feeding ratio
# Q1
plot_map_fc2 +
geom_point(data = filter(full_dat, quarter == 1 & prey_group == "tot_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
facet_grid(species2 ~ year) +
scale_color_viridis(trans = "sqrt", name = "") +
labs(title = "Total feeding ratio", subtitle = "Quarter 1") +
theme(legend.text = element_text(angle = 90))
#> filter: removed 13,798 rows (78%), 3,857 rows remaining
ggsave("figures/supp/data_in_space_q1.pdf", width = 17, height = 17, units = "cm")
# Q4
plot_map_fc2 +
geom_point(data = filter(full_dat, quarter == 4 & prey_group == "tot_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
facet_grid(species2 ~ year) +
scale_color_viridis(trans = "sqrt", name = "") +
labs(title = "Total feeding ratio", subtitle = "Quarter 4") +
theme(legend.text = element_text(angle = 90))
#> filter: removed 15,627 rows (89%), 2,028 rows remaining
ggsave("figures/supp/data_in_space_q4.pdf", width = 17, height = 17, units = "cm")
## Benthic feeding ratio
# Q1
plot_map_fc2 +
geom_point(data = filter(full_dat, quarter == 1 & prey_group == "benthic_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
facet_grid(species2 ~ year) +
scale_color_viridis(trans = "sqrt") +
labs(title = "Benthic feeding ratio", subtitle = "Quarter 1") +
theme(legend.text = element_text(angle = 90))
#> filter: removed 13,798 rows (78%), 3,857 rows remaining
# Q4
plot_map_fc2 +
geom_point(data = filter(full_dat, quarter == 4 & prey_group == "benthic_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
facet_grid(species2 ~ year) +
scale_color_viridis(trans = "sqrt") +
labs(title = "Benthic feeding ratio", subtitle = "Quarter 4") +
theme(legend.text = element_text(angle = 90))
#> filter: removed 15,627 rows (89%), 2,028 rows remaining
## Saduria feeding ratio
# Q1
plot_map_fc2 +
geom_point(data = filter(full_dat, quarter == 1 & prey_group == "saduria_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
facet_grid(species2 ~ year) +
scale_color_viridis(trans = "sqrt") +
labs(title = "Saduria feeding ratio", subtitle = "Quarter 1") +
theme(legend.text = element_text(angle = 90))
#> filter: removed 13,798 rows (78%), 3,857 rows remaining
# Q4
plot_map_fc2 +
geom_point(data = filter(full_dat, quarter == 4 & prey_group == "saduria_feeding_ratio"), aes(X*1000, Y*1000, color = feeding_ratio)) +
facet_grid(species2 ~ year) +
scale_color_viridis(trans = "sqrt") +
labs(title = "Saduria feeding ratio", subtitle = "Quarter 4") +
theme(legend.text = element_text(angle = 90))
#> filter: removed 15,627 rows (89%), 2,028 rows remaining
Plot data (response vs explanatory variables)
ggplot(full_dat, aes(depth, feeding_ratio)) +
geom_point() +
stat_smooth() +
facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
ggplot(full_dat, aes(oxy, feeding_ratio)) +
geom_point() +
stat_smooth() +
facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
ggplot(full_dat, aes(temp, feeding_ratio)) +
geom_point() +
stat_smooth() +
facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
ggplot(full_dat, aes(density_flounder_sc, feeding_ratio)) +
geom_point() +
stat_smooth() +
facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
ggplot(full_dat, aes(density_cod_large_sc, feeding_ratio)) +
geom_point() +
stat_smooth() +
facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
ggplot(full_dat, aes(density_cod_small_sc, feeding_ratio)) +
geom_point() +
stat_smooth() +
facet_wrap(prey_group ~ species2, scales = "free")
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
# Saduria
ggplot(full_dat %>% filter(prey_group == "saduria_feeding_ratio"),
aes(density_saduria_sc, feeding_ratio)) +
geom_point() +
stat_smooth() +
facet_wrap(~ species2, scales = "free", ncol = 2)
#> filter: removed 11,770 rows (67%), 5,885 rows remaining
#> `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Fit spatiotemporal model to feeding ratio
sdmTMB models with densities as covariates to stomach content data. First make mesh:
# Small cod
small_cod_spde <- make_mesh(small_cod_stomach,
c("X", "Y"),
cutoff = 15,
type = "kmeans",
seed = 42)
plot(small_cod_spde)
# Large cod
large_cod_spde <- make_mesh(large_cod_stomach,
c("X", "Y"),
n_knots = 75,
type = "kmeans",
seed = 42)
plot(large_cod_spde)
# Flounder
flounder_spde <- make_mesh(flounder_stomach,
c("X", "Y"),
n_knots = 75,
type = "kmeans",
seed = 42)
plot(flounder_spde)
# s_cod_tot0 <- sdmTMB(
# data = small_cod_stomach,
# formula = tot_feeding_ratio ~ 0 + fyear + quarter,
# family = tweedie(link = "log"),
# time = "year",
# spatiotemporal = "off",
# spatial = "off",
# mesh = small_cod_spde,
# control = sdmTMBcontrol(newton_loops = 1)
# )
#
# tidy(s_cod_tot0)
# sanity(s_cod_tot0)
#
# # Plot residuals in space
# small_cod_stomach$tot_resid0 <- residuals(s_cod_tot0)
#
# plot_map_fc2 +
# geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid0), size = 3,
# shape = 21, color = "grey50", stroke = 0.2) +
# scale_fill_gradient2()
#
#
# s_cod_tot <- sdmTMB(
# data = small_cod_stomach,
# formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
# family = tweedie(link = "log"),
# time = "year",
# spatiotemporal = "off",
# spatial = "off",
# mesh = small_cod_spde,
# priors = sdmTMBpriors(
# b = normal(c(rep(NA, 7), rep(0, 7)),
# c(rep(NA, 7), rep(1, 7)))),
# control = sdmTMBcontrol(newton_loops = 1)
# )
#
# tidy(s_cod_tot)
# sanity(s_cod_tot)
#
# # Plot residuals in space
# small_cod_stomach$tot_resid <- residuals(s_cod_tot)
#
# # plot_map_fc2 +
# # geom_point(data = small_cod_stomach, aes(X*1000, Y*1000, fill = tot_resid), size = 3,
# # shape = 21, color = "grey50", stroke = 0.2) +
# # facet_wrap(~ year) +
# # scale_fill_gradient2()
# #
# # plot_map_fc2 +
# # geom_point(data = filter(small_cod_stomach, tot_resid < 5), aes(X*1000, Y*1000, fill = tot_resid), size = 3,
# # shape = 21, color = "grey50", stroke = 0.2) +
# # facet_wrap(~ year) +
# # scale_fill_gradient2()
#
# plot_map_fc2 +
# geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid), size = 3,
# shape = 21, color = "grey50", stroke = 0.2) +
# scale_fill_gradient2()
#
# # Fit model now with spatial random effect
# s_cod_tot2 <- sdmTMB(
# data = small_cod_stomach,
# formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
# family = tweedie(link = "log"),
# time = "year",
# spatiotemporal = "off",
# spatial = "on",
# mesh = small_cod_spde,
# priors = sdmTMBpriors(
# b = normal(c(rep(NA, 7), rep(0, 7)),
# c(rep(NA, 7), rep(1, 7)))),
# control = sdmTMBcontrol(newton_loops = 1)
# )
#
# print(s_cod_tot2)
# tidy(s_cod_tot2)
# sanity(s_cod_tot2)
#
# small_cod_stomach$tot_resid2 <- residuals(s_cod_tot2)
#
# plot_map_fc2 +
# geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid2), size = 3,
# shape = 21, color = "grey50", stroke = 0.2) +
# scale_fill_gradient2()
#
# # Fit model now with spatiotemporal random effect
# s_cod_tot3 <- sdmTMB(
# data = small_cod_stomach,
# formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
# family = tweedie(link = "log"),
# time = "year",
# spatiotemporal = "IID",
# spatial = "off",
# mesh = small_cod_spde,
# priors = sdmTMBpriors(
# b = normal(c(rep(NA, 7), rep(0, 7)),
# c(rep(NA, 7), rep(1, 7)))),
# control = sdmTMBcontrol(newton_loops = 1)
# )
#
# print(s_cod_tot3)
# tidy(s_cod_tot3)
# sanity(s_cod_tot3)
#
# small_cod_stomach$tot_resid3 <- residuals(s_cod_tot3)
#
# plot_map_fc2 +
# geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid3), size = 3,
# shape = 21, color = "grey50", stroke = 0.2) +
# scale_fill_gradient2()
#
# # Fit model now with spatiotemporal AND spatial random effect
# s_cod_tot4 <- sdmTMB(
# data = small_cod_stomach,
# formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc,
# family = tweedie(link = "log"),
# time = "year",
# spatiotemporal = "IID",
# spatial = "on",
# mesh = small_cod_spde,
# priors = sdmTMBpriors(
# b = normal(c(rep(NA, 7), rep(0, 7)),
# c(rep(NA, 7), rep(1, 7)))),
# control = sdmTMBcontrol(newton_loops = 1)
# )
#
# print(s_cod_tot4)
# tidy(s_cod_tot4)
# sanity(s_cod_tot4)
#
# small_cod_stomach$tot_resid4 <- residuals(s_cod_tot4)
#
# plot_map_fc2 +
# geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid4), size = 3,
# shape = 21, color = "grey50", stroke = 0.2) +
# scale_fill_gradient2()
#
#
# # Just random effects...
# small_cod_stomach$ices_rect <- as.factor(small_cod_stomach$ices_rect)
#
# s_cod_tot5 <- sdmTMB(
# data = small_cod_stomach,
# formula = tot_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*oxy_sc + (1|ices_rect) + (1|haul_id),
# family = tweedie(link = "log"),
# time = "year",
# spatiotemporal = "off",
# spatial = "off",
# mesh = small_cod_spde,
# priors = sdmTMBpriors(
# b = normal(c(rep(NA, 7), rep(0, 7)),
# c(rep(NA, 7), rep(1, 7)))),
# control = sdmTMBcontrol(newton_loops = 1)
# )
#
# print(s_cod_tot5)
# tidy(s_cod_tot5)
# sanity(s_cod_tot5)
#
# small_cod_stomach$tot_resid5 <- residuals(s_cod_tot5)
#
# plot_map_fc2 +
# geom_point(data = filter(small_cod_stomach, tot_resid < 4), aes(X*1000, Y*1000, fill = tot_resid5), size = 3,
# shape = 21, color = "grey50", stroke = 0.2) +
# scale_fill_gradient2()
#
#
# # Plot together
# long <- small_cod_stomach %>%
# pivot_longer(c(tot_resid0, tot_resid, tot_resid2, tot_resid3, tot_resid4, tot_resid5), names_to = "Model", values_to = "Residuals")
#
# plot_map_fc2 +
# geom_jitter(data = long %>% filter(Residuals < 4), aes(X*1000, Y*1000, fill = Residuals), size = 3,
# shape = 21, color = "grey50", stroke = 0.1) +
# scale_fill_gradient2() +
# facet_wrap(~Model, ncol = 3)
#
# plot_map_fc2 +
# geom_jitter(data = long %>% filter(Residuals < 4), aes(X*1000, Y*1000, fill = Residuals), size = 3,
# shape = 21, color = "grey50", stroke = 0.1) +
# scale_fill_gradient2() +
# facet_grid(year ~ Model)
#
#
# # Plot spatial and spatiotemporal random effect
# pred_grid_small <- pred_grid %>%
# filter(quarter == 4) %>%
# mutate(X = X/1000,
# Y = Y/1000) %>%
# filter(year %in% unique(small_cod_stomach$year)) %>%
# filter(X < max(small_cod_stomach$X)) %>%
# filter(X > min(small_cod_stomach$X)) %>%
# filter(Y < max(small_cod_stomach$Y)) %>%
# filter(Y > min(small_cod_stomach$Y))
#
# pred_grid_small <- pred_grid_small %>%
# mutate(fyear = factor(year),
# fquarter = factor(quarter),
# temp_sc = 0,
# depth_sc = 0,
# density_cod_small_sc = 0,
# density_cod_large_sc = 0,
# density_flounder_sc = 0,
# oxy_sc = 0)
#
# spat_pred <- predict(s_cod_tot2, newdata = pred_grid_small)
#
# plot_map_fc2 +
# geom_raster(data = spat_pred, aes(X*1000, Y*1000, fill = omega_s)) +
# scale_fill_gradient2()
#
# spat_pred2 <- predict(s_cod_tot3, newdata = pred_grid_small)
#
# plot_map_fc2 +
# geom_raster(data = spat_pred2, aes(X*1000, Y*1000, fill = epsilon_st)) +
# scale_fill_gradient2() +
# facet_wrap(~year)
#
# AIC(s_cod_tot0, s_cod_tot, s_cod_tot2, s_cod_tot3, s_cod_tot4, s_cod_tot5)Small cod
# Total
s_cod_tot <- sdmTMB(
data = small_cod_stomach,
formula = tot_feeding_ratio ~
0 +
fyear +
fquarter +
temp_sc +
depth_sc +
oxy_sc +
density_cod_small_sc +
density_cod_large_sc +
density_flounder_sc +
(1|haul_id),
family = delta_lognormal(),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = small_cod_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 6)),
c(rep(NA, 7), rep(1, 6)))),
control = sdmTMBcontrol(newton_loops = 1)
)
print(s_cod_tot)
#> Model fit by ML ['sdmTMB']
#> Formula: tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc +
#> Formula: oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc +
#> Formula: (1 | haul_id)
#> Mesh: small_cod_spde
#> Time column: year
#> Data: small_cod_stomach
#> Family: delta_lognormal(link1 = 'logit', link2 = 'log')
#>
#> Delta/hurdle model 1: -----------------------------------
#> Family: binomial(link = 'logit')
#> coef.est coef.se
#> fyear2015 1.91 0.64
#> fyear2016 1.53 0.39
#> fyear2017 2.44 0.35
#> fyear2018 1.64 0.30
#> fyear2019 1.47 0.60
#> fyear2020 0.80 0.23
#> fquarter4 0.06 0.52
#> temp_sc -0.03 0.16
#> depth_sc -0.36 0.18
#> oxy_sc 0.26 0.19
#> density_cod_small_sc 0.01 0.30
#> density_cod_large_sc 0.06 0.29
#> density_flounder_sc 0.20 0.13
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 0.77
#>
#>
#> Delta/hurdle model 2: -----------------------------------
#> Family: lognormal(link = 'log')
#> coef.est coef.se
#> fyear2015 -4.56 0.34
#> fyear2016 -4.92 0.23
#> fyear2017 -4.34 0.18
#> fyear2018 -5.15 0.18
#> fyear2019 -4.53 0.32
#> fyear2020 -5.00 0.16
#> fquarter4 -0.14 0.27
#> temp_sc -0.16 0.10
#> depth_sc -0.06 0.11
#> oxy_sc 0.00 0.12
#> density_cod_small_sc -0.16 0.19
#> density_cod_large_sc 0.11 0.17
#> density_flounder_sc -0.03 0.08
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 0.53
#>
#> Dispersion parameter: 1.26
#>
#> ML criterion at convergence: -3016.195
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(s_cod_tot)
#> term estimate std.error
#> 1 fyear2015 1.90743421 0.6411582
#> 2 fyear2016 1.53353316 0.3900081
#> 3 fyear2017 2.43500762 0.3512151
#> 4 fyear2018 1.64370953 0.2964206
#> 5 fyear2019 1.46943783 0.6031040
#> 6 fyear2020 0.80352813 0.2348412
#> 7 fquarter4 0.05884692 0.5173597
#> 8 temp_sc -0.03332226 0.1620087
#> 9 depth_sc -0.36318086 0.1837213
#> 10 oxy_sc 0.26380179 0.1921616
#> 11 density_cod_small_sc 0.01111117 0.3029440
#> 12 density_cod_large_sc 0.06244649 0.2880924
#> 13 density_flounder_sc 0.19657169 0.1297704
sanity(s_cod_tot)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✖ `b_j2` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
# Benthic
s_cod_ben <- sdmTMB(
data = small_cod_stomach,
formula = benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
family = delta_lognormal(),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = small_cod_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 6)),
c(rep(NA, 7), rep(1, 6)))),
control = sdmTMBcontrol(newton_loops = 1)
)
print(s_cod_ben)
#> Model fit by ML ['sdmTMB']
#> Formula: benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc +
#> Formula: oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc +
#> Formula: (1 | haul_id)
#> Mesh: small_cod_spde
#> Time column: year
#> Data: small_cod_stomach
#> Family: delta_lognormal(link1 = 'logit', link2 = 'log')
#>
#> Delta/hurdle model 1: -----------------------------------
#> Family: binomial(link = 'logit')
#> coef.est coef.se
#> fyear2015 1.57 0.64
#> fyear2016 1.36 0.39
#> fyear2017 2.15 0.33
#> fyear2018 1.53 0.29
#> fyear2019 0.91 0.60
#> fyear2020 0.68 0.24
#> fquarter4 0.42 0.51
#> temp_sc -0.06 0.16
#> depth_sc -0.38 0.18
#> oxy_sc 0.39 0.19
#> density_cod_small_sc 0.15 0.30
#> density_cod_large_sc -0.12 0.27
#> density_flounder_sc 0.22 0.13
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 0.78
#>
#>
#> Delta/hurdle model 2: -----------------------------------
#> Family: lognormal(link = 'log')
#> coef.est coef.se
#> fyear2015 -4.83 0.34
#> fyear2016 -5.14 0.23
#> fyear2017 -4.54 0.18
#> fyear2018 -5.23 0.18
#> fyear2019 -4.95 0.32
#> fyear2020 -5.21 0.16
#> fquarter4 0.05 0.27
#> temp_sc -0.15 0.10
#> depth_sc -0.05 0.11
#> oxy_sc 0.15 0.12
#> density_cod_small_sc -0.03 0.18
#> density_cod_large_sc -0.04 0.17
#> density_flounder_sc 0.00 0.08
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 0.52
#>
#> Dispersion parameter: 1.22
#>
#> ML criterion at convergence: -2988.002
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(s_cod_ben)
#> term estimate std.error
#> 1 fyear2015 1.57305072 0.6372890
#> 2 fyear2016 1.35679232 0.3871899
#> 3 fyear2017 2.14593615 0.3327014
#> 4 fyear2018 1.52583561 0.2940587
#> 5 fyear2019 0.91017699 0.5955846
#> 6 fyear2020 0.68246213 0.2358374
#> 7 fquarter4 0.42191915 0.5084646
#> 8 temp_sc -0.05686552 0.1614583
#> 9 depth_sc -0.38373605 0.1808483
#> 10 oxy_sc 0.38603293 0.1892451
#> 11 density_cod_small_sc 0.15307617 0.2985473
#> 12 density_cod_large_sc -0.12246010 0.2670390
#> 13 density_flounder_sc 0.22236500 0.1304145
sanity(s_cod_ben)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j2` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
# Saduria
s_cod_sad <- sdmTMB(
data = small_cod_stomach,
formula = saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + oxy_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*density_saduria_sc + (1|haul_id),
family = delta_lognormal(),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = small_cod_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 8)),
c(rep(NA, 7), rep(1, 8)))),
control = sdmTMBcontrol(newton_loops = 1)
)
print(s_cod_sad)
#> Model fit by ML ['sdmTMB']
#> Formula: saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + oxy_sc +
#> Formula: depth_sc + density_cod_small_sc + density_cod_large_sc +
#> Formula: density_flounder_sc * density_saduria_sc + (1 | haul_id)
#> Mesh: small_cod_spde
#> Time column: year
#> Data: small_cod_stomach
#> Family: delta_lognormal(link1 = 'logit', link2 = 'log')
#>
#> Delta/hurdle model 1: -----------------------------------
#> Family: binomial(link = 'logit')
#> coef.est coef.se
#> fyear2015 -1.97 0.83
#> fyear2016 -2.33 0.55
#> fyear2017 -2.89 0.47
#> fyear2018 -2.75 0.48
#> fyear2019 -3.49 0.87
#> fyear2020 -2.46 0.37
#> fquarter4 -0.20 0.68
#> temp_sc 0.49 0.23
#> oxy_sc 0.57 0.28
#> depth_sc 0.41 0.29
#> density_cod_small_sc 0.11 0.40
#> density_cod_large_sc 0.02 0.39
#> density_flounder_sc -1.12 0.36
#> density_saduria_sc 0.55 0.23
#> density_flounder_sc:density_saduria_sc 0.59 0.43
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 1.14
#>
#>
#> Delta/hurdle model 2: -----------------------------------
#> Family: lognormal(link = 'log')
#> coef.est coef.se
#> fyear2015 -6.19 0.69
#> fyear2016 -5.20 0.47
#> fyear2017 -6.21 0.42
#> fyear2018 -5.17 0.42
#> fyear2019 -5.86 0.78
#> fyear2020 -4.87 0.33
#> fquarter4 0.17 0.57
#> temp_sc -0.20 0.21
#> oxy_sc -0.16 0.25
#> depth_sc 0.07 0.32
#> density_cod_small_sc -0.18 0.50
#> density_cod_large_sc 0.27 0.66
#> density_flounder_sc 0.20 0.36
#> density_saduria_sc 0.27 0.21
#> density_flounder_sc:density_saduria_sc 0.30 0.42
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 0.14
#>
#> Dispersion parameter: 1.30
#>
#> ML criterion at convergence: -271.004
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(s_cod_sad)
#> term estimate std.error
#> 1 fyear2015 -1.96975855 0.8299213
#> 2 fyear2016 -2.33073173 0.5472003
#> 3 fyear2017 -2.88815853 0.4729243
#> 4 fyear2018 -2.75489048 0.4813353
#> 5 fyear2019 -3.49083692 0.8668283
#> 6 fyear2020 -2.45756405 0.3749937
#> 7 fquarter4 -0.19681017 0.6823595
#> 8 temp_sc 0.49499626 0.2292552
#> 9 oxy_sc 0.56663816 0.2773261
#> 10 depth_sc 0.40580797 0.2873256
#> 11 density_cod_small_sc 0.10921056 0.4002075
#> 12 density_cod_large_sc 0.01993056 0.3870146
#> 13 density_flounder_sc -1.11983774 0.3647768
#> 14 density_saduria_sc 0.55029773 0.2316145
#> 15 density_flounder_sc:density_saduria_sc 0.58664443 0.4258442
sanity(s_cod_sad)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100# Test residuals for single model, tweedie vs lognormal
# s_cod_tot2 <- sdmTMB(
# data = small_cod_stomach,
# formula = tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + density_cod_small_sc +
# density_cod_large_sc + density_flounder_sc*oxy_sc + (1|haul_id),
# family = tweedie(),
# time = "year",
# spatiotemporal = "off",
# spatial = "off",
# mesh = small_cod_spde,
# priors = sdmTMBpriors(
# b = normal(c(rep(NA, 7), rep(0, 7)),
# c(rep(NA, 7), rep(1, 7)))),
# control = sdmTMBcontrol(newton_loops = 1)
# )
#
# # Delta-Lognormal
# small_cod_stomach$res <- residuals(s_cod_tot)
#
# p1 <- plot_map_fc2 +
# geom_point(data = small_cod_stomach, aes(X*1000, Y*1000, fill = res), size = 2, shape = 21, color = "grey50", stroke = 0.2) +
# scale_fill_gradient2() +
# ggtitle("Delta-Lognormal")
#
# mcmc_res <- residuals(s_cod_tot, type = "mle-mcmc", mcmc_iter = 201, mcmc_warmup = 200)
# p1b <- ggplot(as.data.frame(mcmc_res), aes(sample = mcmc_res)) + stat_qq() + stat_qq_line() + theme(aspect.ratio = 1)
# qqnorm(mcmc_res);qqline(mcmc_res)
#
# # Tweedie
# small_cod_stomach$res2 <- residuals(s_cod_tot2)
#
# p2 <- plot_map_fc2 +
# geom_point(data = small_cod_stomach, aes(X*1000, Y*1000, fill = res2), size = 2, shape = 21, color = "grey50", stroke = 0.2) +
# scale_fill_gradient2() +
# ggtitle("Tweedie")
#
# mcmc_res2 <- residuals(s_cod_tot2, type = "mle-mcmc", mcmc_iter = 201, mcmc_warmup = 200) %>% as.data.frame()
# p2b <- ggplot(mcmc_res2, aes(sample = V1)) + stat_qq() + stat_qq_line() + theme(aspect.ratio = 1)
#
# # Plot!
# (p1 + p2 + p1b + p2b) + plot_layout(heights = c(2, 1))Large cod
# Total
l_cod_tot <- sdmTMB(
data = large_cod_stomach,
formula = tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
family = delta_lognormal(),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = large_cod_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 6)),
c(rep(NA, 7), rep(1, 6)))),
control = sdmTMBcontrol(newton_loops = 1)
)
print(l_cod_tot)
#> Model fit by ML ['sdmTMB']
#> Formula: tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc +
#> Formula: oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc +
#> Formula: (1 | haul_id)
#> Mesh: large_cod_spde
#> Time column: year
#> Data: large_cod_stomach
#> Family: delta_lognormal(link1 = 'logit', link2 = 'log')
#>
#> Delta/hurdle model 1: -----------------------------------
#> Family: binomial(link = 'logit')
#> coef.est coef.se
#> fyear2015 1.38 0.35
#> fyear2016 1.49 0.23
#> fyear2017 1.96 0.20
#> fyear2018 1.34 0.16
#> fyear2019 1.06 0.35
#> fyear2020 1.22 0.17
#> fquarter4 0.38 0.31
#> temp_sc -0.46 0.11
#> depth_sc -0.11 0.13
#> oxy_sc 0.05 0.13
#> density_cod_small_sc 0.76 0.30
#> density_cod_large_sc -0.60 0.30
#> density_flounder_sc 0.00 0.08
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 0.44
#>
#>
#> Delta/hurdle model 2: -----------------------------------
#> Family: lognormal(link = 'log')
#> coef.est coef.se
#> fyear2015 -3.08 0.33
#> fyear2016 -3.80 0.21
#> fyear2017 -3.44 0.17
#> fyear2018 -3.58 0.16
#> fyear2019 -3.71 0.32
#> fyear2020 -3.74 0.17
#> fquarter4 -0.79 0.27
#> temp_sc 0.04 0.10
#> depth_sc 0.03 0.13
#> oxy_sc -0.05 0.12
#> density_cod_small_sc 0.20 0.28
#> density_cod_large_sc -0.05 0.28
#> density_flounder_sc -0.09 0.08
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 0.55
#>
#> Dispersion parameter: 1.76
#>
#> ML criterion at convergence: -4845.201
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(l_cod_tot)
#> term estimate std.error
#> 1 fyear2015 1.376267080 0.35248671
#> 2 fyear2016 1.489466984 0.23264291
#> 3 fyear2017 1.955719828 0.20279435
#> 4 fyear2018 1.342176159 0.16264333
#> 5 fyear2019 1.063579950 0.35163094
#> 6 fyear2020 1.224150901 0.17030026
#> 7 fquarter4 0.382795458 0.31010761
#> 8 temp_sc -0.459569859 0.11054982
#> 9 depth_sc -0.109393256 0.13361705
#> 10 oxy_sc 0.054455940 0.12583284
#> 11 density_cod_small_sc 0.757360174 0.30113449
#> 12 density_cod_large_sc -0.597163459 0.29919650
#> 13 density_flounder_sc 0.002566405 0.07747939
sanity(l_cod_tot)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
# Benthic
l_cod_ben <- sdmTMB(
data = large_cod_stomach,
formula = benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
family = delta_lognormal(),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = large_cod_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 6)),
c(rep(NA, 7), rep(1, 6)))),
control = sdmTMBcontrol(newton_loops = 1)
)
print(l_cod_ben)
#> Model fit by ML ['sdmTMB']
#> Formula: benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc +
#> Formula: oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc +
#> Formula: (1 | haul_id)
#> Mesh: large_cod_spde
#> Time column: year
#> Data: large_cod_stomach
#> Family: delta_lognormal(link1 = 'logit', link2 = 'log')
#>
#> Delta/hurdle model 1: -----------------------------------
#> Family: binomial(link = 'logit')
#> coef.est coef.se
#> fyear2015 -0.18 0.36
#> fyear2016 0.15 0.22
#> fyear2017 0.41 0.18
#> fyear2018 -0.09 0.16
#> fyear2019 -0.78 0.34
#> fyear2020 0.02 0.17
#> fquarter4 1.26 0.30
#> temp_sc -0.38 0.11
#> depth_sc -0.10 0.14
#> oxy_sc 0.28 0.13
#> density_cod_small_sc 0.25 0.30
#> density_cod_large_sc -0.21 0.29
#> density_flounder_sc 0.05 0.08
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 0.58
#>
#>
#> Delta/hurdle model 2: -----------------------------------
#> Family: lognormal(link = 'log')
#> coef.est coef.se
#> fyear2015 -5.50 0.31
#> fyear2016 -5.71 0.20
#> fyear2017 -5.56 0.17
#> fyear2018 -5.63 0.16
#> fyear2019 -6.35 0.30
#> fyear2020 -5.72 0.17
#> fquarter4 0.39 0.26
#> temp_sc -0.12 0.10
#> depth_sc -0.04 0.13
#> oxy_sc 0.08 0.12
#> density_cod_small_sc -0.41 0.27
#> density_cod_large_sc 0.52 0.27
#> density_flounder_sc -0.09 0.07
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 0.49
#>
#> Dispersion parameter: 1.55
#>
#> ML criterion at convergence: -4771.129
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(l_cod_ben)
#> term estimate std.error
#> 1 fyear2015 -0.17736775 0.35712898
#> 2 fyear2016 0.15284441 0.21906882
#> 3 fyear2017 0.41442854 0.18416163
#> 4 fyear2018 -0.09346653 0.15820153
#> 5 fyear2019 -0.78207442 0.34474478
#> 6 fyear2020 0.01629230 0.17278533
#> 7 fquarter4 1.26366931 0.29920972
#> 8 temp_sc -0.38340307 0.10941454
#> 9 depth_sc -0.09719543 0.13672989
#> 10 oxy_sc 0.28408572 0.12649178
#> 11 density_cod_small_sc 0.25424771 0.29604417
#> 12 density_cod_large_sc -0.20528910 0.29446028
#> 13 density_flounder_sc 0.04657151 0.08037281
sanity(l_cod_ben)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
# Saduria
l_cod_sad <- sdmTMB(
data = large_cod_stomach,
formula = saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + oxy_sc + depth_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc*density_saduria_sc + (1|haul_id),
family = delta_lognormal(),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = large_cod_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 8)),
c(rep(NA, 7), rep(1, 8)))),
control = sdmTMBcontrol(newton_loops = 1)
)
print(l_cod_sad)
#> Model fit by ML ['sdmTMB']
#> Formula: saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + oxy_sc +
#> Formula: depth_sc + density_cod_small_sc + density_cod_large_sc +
#> Formula: density_flounder_sc * density_saduria_sc + (1 | haul_id)
#> Mesh: large_cod_spde
#> Time column: year
#> Data: large_cod_stomach
#> Family: delta_lognormal(link1 = 'logit', link2 = 'log')
#>
#> Delta/hurdle model 1: -----------------------------------
#> Family: binomial(link = 'logit')
#> coef.est coef.se
#> fyear2015 -1.80 0.83
#> fyear2016 -3.39 0.54
#> fyear2017 -3.64 0.51
#> fyear2018 -3.65 0.44
#> fyear2019 -4.15 0.92
#> fyear2020 -2.82 0.43
#> fquarter4 -0.05 0.71
#> temp_sc 0.00 0.25
#> oxy_sc 0.39 0.31
#> depth_sc 0.23 0.36
#> density_cod_small_sc 0.33 0.56
#> density_cod_large_sc -0.38 0.60
#> density_flounder_sc -0.83 0.31
#> density_saduria_sc 0.40 0.23
#> density_flounder_sc:density_saduria_sc 0.40 0.38
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 1.53
#>
#>
#> Delta/hurdle model 2: -----------------------------------
#> Family: lognormal(link = 'log')
#> coef.est coef.se
#> fyear2015 -5.58 0.54
#> fyear2016 -6.02 0.35
#> fyear2017 -5.74 0.36
#> fyear2018 -5.95 0.29
#> fyear2019 -6.02 0.71
#> fyear2020 -6.07 0.32
#> fquarter4 -0.78 0.53
#> temp_sc 0.53 0.18
#> oxy_sc 0.33 0.24
#> depth_sc 0.35 0.31
#> density_cod_small_sc -0.16 0.50
#> density_cod_large_sc 0.28 0.77
#> density_flounder_sc -0.12 0.22
#> density_saduria_sc 0.16 0.18
#> density_flounder_sc:density_saduria_sc -0.12 0.34
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 0.39
#>
#> Dispersion parameter: 1.29
#>
#> ML criterion at convergence: -498.453
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(l_cod_sad)
#> term estimate std.error
#> 1 fyear2015 -1.799469107 0.8330093
#> 2 fyear2016 -3.385412445 0.5368501
#> 3 fyear2017 -3.640778227 0.5114008
#> 4 fyear2018 -3.650167817 0.4403155
#> 5 fyear2019 -4.147620658 0.9194658
#> 6 fyear2020 -2.821250957 0.4257319
#> 7 fquarter4 -0.052218529 0.7129595
#> 8 temp_sc 0.002607801 0.2518307
#> 9 oxy_sc 0.389346476 0.3110741
#> 10 depth_sc 0.227847260 0.3615738
#> 11 density_cod_small_sc 0.328026403 0.5551598
#> 12 density_cod_large_sc -0.382032051 0.5989242
#> 13 density_flounder_sc -0.834920812 0.3136644
#> 14 density_saduria_sc 0.397871281 0.2278797
#> 15 density_flounder_sc:density_saduria_sc 0.396319690 0.3848843
sanity(l_cod_sad)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100Flounder
# Total
fle_tot <- sdmTMB(
data = flounder_stomach,
formula = tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
family = delta_lognormal(),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = flounder_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 6)),
c(rep(NA, 7), rep(1, 6)))),
control = sdmTMBcontrol(newton_loops = 1)
)
print(fle_tot)
#> Model fit by ML ['sdmTMB']
#> Formula: tot_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc +
#> Formula: oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc +
#> Formula: (1 | haul_id)
#> Mesh: flounder_spde
#> Time column: year
#> Data: flounder_stomach
#> Family: delta_lognormal(link1 = 'logit', link2 = 'log')
#>
#> Delta/hurdle model 1: -----------------------------------
#> Family: binomial(link = 'logit')
#> coef.est coef.se
#> fyear2015 0.53 0.53
#> fyear2016 0.30 0.39
#> fyear2017 0.36 0.25
#> fyear2018 0.12 0.25
#> fyear2019 -0.03 0.47
#> fyear2020 1.42 0.23
#> fquarter4 0.86 0.38
#> temp_sc -0.10 0.14
#> depth_sc -0.20 0.17
#> oxy_sc 0.33 0.19
#> density_cod_small_sc 0.63 0.37
#> density_cod_large_sc -0.58 0.37
#> density_flounder_sc 0.03 0.12
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 1.04
#>
#>
#> Delta/hurdle model 2: -----------------------------------
#> Family: lognormal(link = 'log')
#> coef.est coef.se
#> fyear2015 -4.37 0.24
#> fyear2016 -4.40 0.18
#> fyear2017 -3.98 0.13
#> fyear2018 -4.04 0.13
#> fyear2019 -4.06 0.22
#> fyear2020 -3.64 0.10
#> fquarter4 0.37 0.18
#> temp_sc 0.18 0.07
#> depth_sc 0.00 0.09
#> oxy_sc 0.13 0.10
#> density_cod_small_sc -0.19 0.19
#> density_cod_large_sc 0.21 0.19
#> density_flounder_sc -0.07 0.05
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 0.38
#>
#> Dispersion parameter: 1.26
#>
#> ML criterion at convergence: -3530.324
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(fle_tot)
#> term estimate std.error
#> 1 fyear2015 0.52611178 0.5277953
#> 2 fyear2016 0.29633519 0.3855055
#> 3 fyear2017 0.36424865 0.2499490
#> 4 fyear2018 0.11872342 0.2529731
#> 5 fyear2019 -0.03306316 0.4684866
#> 6 fyear2020 1.42321774 0.2260591
#> 7 fquarter4 0.86332478 0.3823155
#> 8 temp_sc -0.10027050 0.1439928
#> 9 depth_sc -0.20081704 0.1677329
#> 10 oxy_sc 0.32846294 0.1869042
#> 11 density_cod_small_sc 0.62940403 0.3726341
#> 12 density_cod_large_sc -0.57761745 0.3693939
#> 13 density_flounder_sc 0.03169148 0.1150707
sanity(fle_tot)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✖ `b_j2` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
# Benthic
fle_ben <- sdmTMB(
data = flounder_stomach,
formula = benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc + oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
family = delta_lognormal(),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = flounder_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 6)),
c(rep(NA, 7), rep(1, 6)))),
control = sdmTMBcontrol(newton_loops = 1)
)
print(fle_ben)
#> Model fit by ML ['sdmTMB']
#> Formula: benthic_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + depth_sc +
#> Formula: oxy_sc + density_cod_small_sc + density_cod_large_sc + density_flounder_sc +
#> Formula: (1 | haul_id)
#> Mesh: flounder_spde
#> Time column: year
#> Data: flounder_stomach
#> Family: delta_lognormal(link1 = 'logit', link2 = 'log')
#>
#> Delta/hurdle model 1: -----------------------------------
#> Family: binomial(link = 'logit')
#> coef.est coef.se
#> fyear2015 0.51 0.54
#> fyear2016 0.27 0.39
#> fyear2017 0.32 0.25
#> fyear2018 0.14 0.26
#> fyear2019 -0.28 0.48
#> fyear2020 1.33 0.23
#> fquarter4 0.86 0.39
#> temp_sc -0.05 0.15
#> depth_sc -0.22 0.17
#> oxy_sc 0.35 0.19
#> density_cod_small_sc 0.66 0.38
#> density_cod_large_sc -0.62 0.37
#> density_flounder_sc 0.05 0.12
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 1.06
#>
#>
#> Delta/hurdle model 2: -----------------------------------
#> Family: lognormal(link = 'log')
#> coef.est coef.se
#> fyear2015 -4.39 0.23
#> fyear2016 -4.43 0.17
#> fyear2017 -3.98 0.13
#> fyear2018 -4.01 0.13
#> fyear2019 -4.14 0.22
#> fyear2020 -3.65 0.10
#> fquarter4 0.37 0.18
#> temp_sc 0.19 0.07
#> depth_sc -0.07 0.09
#> oxy_sc 0.09 0.10
#> density_cod_small_sc -0.21 0.18
#> density_cod_large_sc 0.23 0.18
#> density_flounder_sc -0.05 0.05
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 0.35
#>
#> Dispersion parameter: 1.28
#>
#> ML criterion at convergence: -3476.879
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(fle_ben)
#> term estimate std.error
#> 1 fyear2015 0.50568356 0.5372639
#> 2 fyear2016 0.27473734 0.3927066
#> 3 fyear2017 0.31836479 0.2545055
#> 4 fyear2018 0.13679684 0.2572819
#> 5 fyear2019 -0.27996769 0.4768295
#> 6 fyear2020 1.32914073 0.2284407
#> 7 fquarter4 0.86119225 0.3891565
#> 8 temp_sc -0.04676748 0.1462127
#> 9 depth_sc -0.21504409 0.1708393
#> 10 oxy_sc 0.35217170 0.1903804
#> 11 density_cod_small_sc 0.65658216 0.3776414
#> 12 density_cod_large_sc -0.61939795 0.3746410
#> 13 density_flounder_sc 0.05410617 0.1171719
sanity(fle_ben)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
# Saduria
fle_sad <- sdmTMB(
data = flounder_stomach,
formula = saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + oxy_sc + depth_sc + density_flounder_sc + density_cod_large_sc + density_cod_small_sc*density_saduria_sc + (1|haul_id),
family = delta_lognormal(),
time = "year",
spatiotemporal = "off",
spatial = "off",
mesh = flounder_spde,
priors = sdmTMBpriors(
b = normal(c(rep(NA, 7), rep(0, 8)),
c(rep(NA, 7), rep(1, 8)))),
control = sdmTMBcontrol(newton_loops = 1)
)
print(fle_sad)
#> Model fit by ML ['sdmTMB']
#> Formula: saduria_feeding_ratio ~ 0 + fyear + fquarter + temp_sc + oxy_sc +
#> Formula: depth_sc + density_flounder_sc + density_cod_large_sc + density_cod_small_sc *
#> Formula: density_saduria_sc + (1 | haul_id)
#> Mesh: flounder_spde
#> Time column: year
#> Data: flounder_stomach
#> Family: delta_lognormal(link1 = 'logit', link2 = 'log')
#>
#> Delta/hurdle model 1: -----------------------------------
#> Family: binomial(link = 'logit')
#> coef.est coef.se
#> fyear2015 -1.08 0.85
#> fyear2016 -0.50 0.65
#> fyear2017 -2.43 0.45
#> fyear2018 -3.65 0.49
#> fyear2019 -2.69 0.82
#> fyear2020 -2.07 0.40
#> fquarter4 -0.31 0.64
#> temp_sc -0.20 0.23
#> oxy_sc 0.26 0.30
#> depth_sc 0.07 0.27
#> density_flounder_sc -0.42 0.26
#> density_cod_large_sc 0.57 0.54
#> density_cod_small_sc -0.82 0.58
#> density_saduria_sc 0.68 0.26
#> density_cod_small_sc:density_saduria_sc -0.11 0.52
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 1.75
#>
#>
#> Delta/hurdle model 2: -----------------------------------
#> Family: lognormal(link = 'log')
#> coef.est coef.se
#> fyear2015 -3.92 0.41
#> fyear2016 -4.00 0.30
#> fyear2017 -5.18 0.24
#> fyear2018 -4.79 0.29
#> fyear2019 -4.09 0.43
#> fyear2020 -4.31 0.22
#> fquarter4 -1.06 0.31
#> temp_sc -0.15 0.13
#> oxy_sc -0.57 0.19
#> depth_sc -0.49 0.20
#> density_flounder_sc -0.39 0.15
#> density_cod_large_sc 0.24 0.36
#> density_cod_small_sc -0.21 0.43
#> density_saduria_sc 0.29 0.15
#> density_cod_small_sc:density_saduria_sc 0.16 0.33
#>
#> Random intercepts:
#> Std. Dev.
#> haul_id 0.53
#>
#> Dispersion parameter: 1.27
#>
#> ML criterion at convergence: -1139.288
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
tidy(fle_sad)
#> term estimate std.error
#> 1 fyear2015 -1.08052722 0.8501067
#> 2 fyear2016 -0.50181784 0.6485335
#> 3 fyear2017 -2.43386977 0.4479369
#> 4 fyear2018 -3.64549688 0.4871173
#> 5 fyear2019 -2.68587056 0.8176042
#> 6 fyear2020 -2.07202390 0.3970650
#> 7 fquarter4 -0.31306727 0.6417368
#> 8 temp_sc -0.19566564 0.2337761
#> 9 oxy_sc 0.26092780 0.3040494
#> 10 depth_sc 0.06517983 0.2721824
#> 11 density_flounder_sc -0.41502450 0.2621885
#> 12 density_cod_large_sc 0.56748090 0.5371357
#> 13 density_cod_small_sc -0.81664754 0.5787022
#> 14 density_saduria_sc 0.68338619 0.2593120
#> 15 density_cod_small_sc:density_saduria_sc -0.10661946 0.5234617
sanity(fle_sad)
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigen values detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No fixed-effect standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100mcmc_res_fle_sad <- residuals(fle_sad,
type = "mle-mcmc",
mcmc_iter = 201,
mcmc_warmup = 200)
#> Constructing atomic log_dbinom_robust
#> Constructing atomic invpd
#>
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.001647 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 16.47 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
#> Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
#> Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
#> Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
#> Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 14.2302 seconds (Warm-up)
#> Chain 1: 0.026116 seconds (Sampling)
#> Chain 1: 14.2563 seconds (Total)
#> Chain 1:
qqnorm(mcmc_res_fle_sad); qqline(mcmc_res_fle_sad)
Plot coefficients
#### Total prey
# small cod total prey
small_cod_tot_pars <- bind_rows(tidy(s_cod_tot, effects = "fixed", conf.int = TRUE, model = 1) %>%
mutate(model = "Binomial"),
tidy(s_cod_tot, effects = "fixed", conf.int = TRUE, model = 2) %>%
mutate(model = "Lognormal")) %>%
mutate(species = "Small cod",
response = "Total prey")
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
# large cod total prey
large_cod_tot_pars <- bind_rows(tidy(l_cod_tot, effects = "fixed", conf.int = TRUE, model = 1) %>%
mutate(model = "Binomial"),
tidy(l_cod_tot, effects = "fixed", conf.int = TRUE, model = 2) %>%
mutate(model = "Lognormal")) %>%
mutate(species = "Large cod",
response = "Total prey")
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
# flounder total prey
flounder_tot_pars <- bind_rows(tidy(fle_tot, effects = "fixed", conf.int = TRUE, model = 1) %>%
mutate(model = "Binomial"),
tidy(fle_tot, effects = "fixed", conf.int = TRUE, model = 2) %>%
mutate(model = "Lognormal")) %>%
mutate(species = "Flounder",
response = "Total prey")
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
#### Benthic prey
# small cod benthic prey
small_cod_ben_pars <- bind_rows(tidy(s_cod_ben, effects = "fixed", conf.int = TRUE, model = 1) %>%
mutate(model = "Binomial"),
tidy(s_cod_ben, effects = "fixed", conf.int = TRUE, model = 2) %>%
mutate(model = "Lognormal")) %>%
mutate(species = "Small cod",
response = "Benthic prey")
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
# large cod benthic prey
large_cod_ben_pars <- bind_rows(tidy(l_cod_ben, effects = "fixed", conf.int = TRUE, model = 1) %>%
mutate(model = "Binomial"),
tidy(l_cod_ben, effects = "fixed", conf.int = TRUE, model = 2) %>%
mutate(model = "Lognormal")) %>%
mutate(species = "Large cod",
response = "Benthic prey")
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
# flounder benthic prey
flounder_ben_pars <- bind_rows(tidy(fle_ben, effects = "fixed", conf.int = TRUE, model = 1) %>%
mutate(model = "Binomial"),
tidy(fle_ben, effects = "fixed", conf.int = TRUE, model = 2) %>%
mutate(model = "Lognormal")) %>%
mutate(species = "Flounder",
response = "Benthic prey")
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
#### Saduria
# small cod saduria
small_cod_sad_pars <- bind_rows(tidy(s_cod_sad, effects = "fixed", conf.int = TRUE, model = 1) %>%
mutate(model = "Binomial"),
tidy(s_cod_sad, effects = "fixed", conf.int = TRUE, model = 2) %>%
mutate(model = "Lognormal")) %>%
mutate(species = "Small cod",
response = "Saduria")
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
# large cod saduria
large_cod_sad_pars <- bind_rows(tidy(l_cod_sad, effects = "fixed", conf.int = TRUE, model = 1) %>%
mutate(model = "Binomial"),
tidy(l_cod_sad, effects = "fixed", conf.int = TRUE, model = 2) %>%
mutate(model = "Lognormal")) %>%
mutate(species = "Large cod",
response = "Saduria")
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
# flounder saduria
flounder_sad_pars <- bind_rows(tidy(fle_sad, effects = "fixed", conf.int = TRUE, model = 1) %>%
mutate(model = "Binomial"),
tidy(fle_sad, effects = "fixed", conf.int = TRUE, model = 2) %>%
mutate(model = "Lognormal")) %>%
mutate(species = "Flounder",
response = "Saduria")
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> new variable 'response' (character) with one unique value and 0% NA
coefs <- bind_rows(small_cod_tot_pars,
large_cod_tot_pars,
flounder_tot_pars,
small_cod_ben_pars,
large_cod_ben_pars,
flounder_ben_pars,
small_cod_sad_pars,
large_cod_sad_pars,
flounder_sad_pars)
coefs <- coefs %>%
mutate(term2 = term,
term2 = ifelse(term == "temp_sc", "Temperature", term2),
term2 = ifelse(term == "oxy_sc", "Oxygen", term2),
term2 = ifelse(term == "depth_sc", "Depth", term2),
term2 = ifelse(term == "density_saduria_sc", "Saduria", term2),
term2 = ifelse(term == "density_cod_large_sc", "Cod (> 29)", term2),
term2 = ifelse(term == "density_cod_small_sc", "Cod (< 29)", term2),
term2 = ifelse(term == "density_cod_small_sc:density_saduria_sc", "Cod (> 29) x Saduria", term2),
term2 = ifelse(term == "density_flounder_sc:density_saduria_sc", "Flounder x Saduria", term2),
term2 = ifelse(term == "density_cod_small_sc:density_saduria_sc", "Small cod x Saduria", term2),
term2 = ifelse(term == "density_flounder_sc", "Flounder", term2))
#> mutate: new variable 'term2' (character) with 16 unique values and 0% NA
coefs <- coefs %>%
mutate(sig = ifelse(estimate > 0 & conf.low > 0, "Y", "N"),
sig = ifelse(estimate < 0 & conf.high < 0, "Y", sig))
#> mutate: new variable 'sig' (character) with 2 unique values and 0% NA
# Fixed continuous effects:
# coefs %>%
# filter(!grepl('year', term)) %>%
# filter(!grepl('quarter', term)) %>%
# ggplot(aes(term2, estimate, color = species)) +
# geom_hline(yintercept = 0, linetype = 2, color = "gray50", size = 0.2) +
# geom_point(position = position_dodge(width = 0.4)) +
# geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0,
# position = position_dodge(width = 0.4), alpha = 0.5) +
# scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
# facet_grid(model~response) +
# labs(x = "Predictor", y = "Standardized coefficient") +
# theme(legend.position = "bottom",
# axis.text.x = element_text(angle = 90),
# aspect.ratio = 1) +
# NULL
#ggsave("figures/fr_coefs.pdf", width = 17, height = 17, units = "cm")
coefs %>%
filter(!grepl('year', term)) %>%
filter(!grepl('quarter', term)) %>%
ggplot(aes(term2, estimate, color = model, fill = model, shape = sig)) +
geom_hline(yintercept = 0, linetype = 2, color = "gray50", size = 0.2) +
geom_point(position = position_dodge(width = 0.4)) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0,
position = position_dodge(width = 0.4), alpha = 0.5) +
scale_color_brewer(palette = "Set1") +
scale_fill_brewer(palette = "Set1") +
scale_shape_manual(values = c(1, 21)) +
guides(fill = "none") +
facet_grid(response ~ species) +
labs(x = "", y = "Standardized coefficient", color = "Model", shape = "95% CI crossing 0") +
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 90, vjust = 0.4)) +
NULL
#> filter: removed 108 rows (44%), 138 rows remaining
#> filter: removed 18 rows (13%), 120 rows remaining
ggsave("figures/fr_coefs.pdf", width = 17, height = 17, units = "cm")
# Year effects
coefs %>%
filter(term %in% c("fyear2015", "fyear2016", "fyear2017", "fyear2018", "fyear2019", "fyear2020", "quarter")) %>%
mutate(term2 = ifelse(term == "fyear2015", "2015", term2),
term2 = ifelse(term == "fyear2016", "2016", term2),
term2 = ifelse(term == "fyear2017", "2017", term2),
term2 = ifelse(term == "fyear2018", "2018", term2),
term2 = ifelse(term == "fyear2019", "2019", term2),
term2 = ifelse(term == "fyear2020", "2020", term2),
term2 = ifelse(term == "quarter", "Quarter", term2)) %>%
filter(!term2 == "Quarter") %>%
ggplot(aes(term2, estimate, color = species)) +
geom_point(position = position_dodge(width = 0.5)) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0,
position = position_dodge(width = 0.5)) +
scale_color_brewer(palette = "Paired", name = "Species", direction = -1) +
facet_grid(model~response) +
labs(x = "Predictor", y = "Standardized coefficient") +
theme(legend.position = "bottom",
axis.text.x = element_text(angle = 90),
aspect.ratio = 1.5) +
NULL
#> filter: removed 138 rows (56%), 108 rows remaining
#> mutate: changed 108 values (100%) of 'term2' (0 new NA)
#> filter: no rows removed
ggsave("figures/supp/fr_coefs_yr.pdf", width = 17, height = 17, units = "cm")# fle_sad <- sdmTMB(
# data = flounder_stomach,
# formula = saduria_feeding_ratio ~ 0 + fyear + quarter + temp_sc + depth_sc + density_cod_small_sc*oxy_sc + density_cod_large_sc + density_flounder_sc + (1|haul_id),
# family = tweedie(link = "log"),
# time = "year",
# spatiotemporal = "off",
# spatial = "off",
# mesh = flounder_spde,
# priors = sdmTMBpriors(
# b = normal(c(rep(NA, 7), rep(0, 7)),
# c(rep(NA, 7), rep(1, 7)))),
# control = sdmTMBcontrol(newton_loops = 1)
# )
# visreg_delta(s_cod_sad, xvar = "density_flounder_sc", model = 1, gg = TRUE, by = "oxy_sc")
# visreg_delta(s_cod_sad, xvar = "density_flounder_sc", model = 2, gg = TRUE)
# visreg_delta(s_cod_sad, xvar = "density_flounder_sc", model = 1, gg = TRUE, by = "oxy_sc")
hist(fle_sad$data$density_cod_small_sc)
fle1 <- visreg_delta(fle_sad,
xvar = "density_cod_small_sc",
model = 1,
plot = FALSE,
by = "density_saduria_sc",
breaks = c(-0.5, 0.5)
)
fle2 <- visreg_delta(fle_sad,
xvar = "density_cod_small_sc",
model = 2,
plot = FALSE,
by = "density_saduria_sc",
breaks = c(-0.5, 0.5))
scod1 <- visreg_delta(s_cod_sad,
xvar = "density_flounder_sc",
model = 1,
plot = FALSE,
by = "density_saduria_sc",
breaks = c(-0.5, 0.5))
scod2 <- visreg_delta(s_cod_sad,
xvar = "density_flounder_sc",
model = 2,
plot = FALSE,
by = "density_saduria_sc",
breaks = c(-0.5, 0.5))
lcod1 <- visreg_delta(l_cod_sad,
xvar = "density_flounder_sc",
model = 1,
plot = FALSE,
by = "density_saduria_sc",
breaks = c(-0.5, 0.5))
lcod2 <- visreg_delta(l_cod_sad,
xvar = "density_flounder_sc",
model = 2,
plot = FALSE,
by = "density_saduria_sc",
breaks = c(-0.5, 0.5))
visreg_fit <- bind_rows(fle1$fit %>% mutate(model = "fle1"),
fle2$fit %>% mutate(model = "fle2"),
scod1$fit %>% mutate(model = "scod1"),
scod2$fit %>% mutate(model = "scod2"),
lcod1$fit %>% mutate(model = "lcod1"),
lcod2$fit %>% mutate(model = "lcod2")
)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
visreg_res <- bind_rows(fle1$res %>% mutate(model = "fle1"),
fle2$res %>% mutate(model = "fle2"),
scod1$res %>% mutate(model = "scod1"),
scod2$res %>% mutate(model = "scod2"),
lcod1$res %>% mutate(model = "lcod1"),
lcod2$res %>% mutate(model = "lcod2")
)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
visreg_fit <- visreg_fit %>%
mutate(species = ifelse(model %in% c("fle1", "fle2"), "Flounder", NA),
species = ifelse(model %in% c("scod1", "scod2"), "Small cod", species),
species = ifelse(model %in% c("lcod1", "lcod2"), "Large cod", species)) %>%
mutate(model_comp = ifelse(model %in% c("fle1", "scod1", "lcod1"), "Binomial", "LogNormal"))
#> mutate: new variable 'species' (character) with 3 unique values and 0% NA
#> mutate: new variable 'model_comp' (character) with 2 unique values and 0% NA
visreg_res <- visreg_res %>%
mutate(species = ifelse(model %in% c("fle1", "fle2"), "Flounder", NA),
species = ifelse(model %in% c("scod1", "scod2"), "Small cod", species),
species = ifelse(model %in% c("lcod1", "lcod2"), "Large cod", species)) %>%
mutate(model_comp = ifelse(model %in% c("fle1", "scod1", "lcod1"), "Binomial", "LogNormal"))
#> mutate: new variable 'species' (character) with 3 unique values and 0% NA
#> mutate: new variable 'model_comp' (character) with 2 unique values and 0% NA
pal <- brewer.pal(n = 9, name = "Set1")[c(2, 9)]
# Flounder
p1 <- ggplot(visreg_fit %>% filter(species == "Flounder"),
aes(x = density_cod_small_sc, y = visregFit, fill = factor(density_saduria_sc), color = factor(density_saduria_sc))) +
facet_grid(model_comp ~ species) +
geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.2, color = NA) +
geom_line() +
labs(x = "Small cod density (scaled)",
color = "Saduria density (scaled)") +
guides(fill = "none") +
geom_point(aes(y = visregRes), data = visreg_res %>% filter(species == "Flounder"), size = 1, alpha = 0.5) +
# scale_color_brewer(palette = "Set1") +
# scale_fill_brewer(palette = "Set1") +
# scale_color_brewer(palette = "Accent") +
# scale_fill_brewer(palette = "Accent") +
scale_color_manual(values = pal) +
scale_fill_manual(values = pal) +
theme(aspect.ratio = 1.1) +
NULL
#> filter: removed 808 rows (67%), 404 rows remaining
#> filter: removed 3,658 rows (54%), 3,174 rows remaining
# Small cod
p2 <- ggplot(visreg_fit %>% filter(species == "Small cod"),
aes(x = density_flounder_sc, y = visregFit, fill = factor(density_saduria_sc), color = factor(density_saduria_sc))) +
#ggh4x::facet_grid2(model_comp ~ species, scales = "free_y", independent = "y") +
facet_grid(model_comp ~ species) +
geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.2, color = NA) +
geom_line() +
labs(x = "Flounder density (scaled)",
color = "Saduria density (scaled)") +
guides(fill = "none") +
geom_point(aes(y = visregRes), data = visreg_res %>% filter(species == "Small cod"), size = 1, alpha = 0.2) +
scale_color_manual(values = pal) +
scale_fill_manual(values = pal) +
theme(aspect.ratio = 1.1) +
NULL
#> filter: removed 808 rows (67%), 404 rows remaining
#> filter: removed 5,536 rows (81%), 1,296 rows remaining
# Large cod
p3 <- ggplot(visreg_fit %>% filter(species == "Large cod"),
aes(x = density_flounder_sc, y = visregFit, fill = factor(density_saduria_sc), color = factor(density_saduria_sc))) +
#ggh4x::facet_grid2(model_comp ~ species, scales = "free_y", independent = "y") +
#ggh4x::facet_grid2(model_comp ~ species, scales = "free_y", independent = "y") +
facet_grid(model_comp ~ species) +
geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.2, color = NA) +
geom_line() +
labs(x = "Flounder density (scaled)",
color = "Saduria density (scaled)") +
guides(fill = "none") +
geom_point(aes(y = visregRes), data = visreg_res %>% filter(species == "Large cod"), size = 1, alpha = 0.2) +
scale_color_manual(values = pal) +
scale_fill_manual(values = pal) +
theme(aspect.ratio = 1.1) +
NULL
#> filter: removed 808 rows (67%), 404 rows remaining
#> filter: removed 4,470 rows (65%), 2,362 rows remaining
(p1 | p2 | p3) + plot_layout(guides = "collect")
# visreg_fit %>%
# #filter(species == "Large cod") %>%
# group_by(species, model) %>%
# summarise(max_covar = max(density_flounder_sc),
# min_covar = min(density_flounder_sc))
# Option 2
# ggplot(visreg_fit, aes(x = density_flounder_sc, y = visregFit, fill = factor(density_saduria_sc), color = factor(density_saduria_sc))) +
# facet_grid2(model_comp ~ species) +
# geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.2, color = NA) +
# geom_line() +
# labs(x = "Flounder density (scaled)",
# color = "Saduria density") +
# guides(fill = "none") +
# geom_point(aes(y = visregRes), data = visreg_res, size = 1, alpha = 0.2) +
# scale_color_brewer(palette ="Dark2") +
# scale_fill_brewer(palette = "Dark2") +
# theme(aspect.ratio = 1.1) +
# NULL
ggsave("figures/conditional_saduria.pdf", width = 17, height = 17, units = "cm")